''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Cephes Math Library Release 2.8:  June, 2000
'Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
'
'Contributors:
'    * Sergey Bochkanov (ALGLIB project). Translation from C to
'      pseudocode.
'
'See subroutines comments for additional copyrights.
'
'>>> SOURCE LICENSE >>>
'This program is free software; you can redistribute it and/or modify
'it under the terms of the GNU General Public License as published by
'the Free Software Foundation (www.fsf.org); either version 2 of the 
'License, or (at your option) any later version.
'
'This program is distributed in the hope that it will be useful,
'but WITHOUT ANY WARRANTY; without even the implied warranty of
'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
'GNU General Public License for more details.
'
'A copy of the GNU General Public License is available at
'http://www.fsf.org/licensing/licenses
'
'>>> END OF LICENSE >>>
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'F distribution
'
'Returns the area from zero to x under the F density
'function (also known as Snedcor's density or the
'variance ratio density).  This is the density
'of x = (u1/df1)/(u2/df2), where u1 and u2 are random
'variables having Chi square distributions with df1
'and df2 degrees of freedom, respectively.
'The incomplete beta integral is used, according to the
'formula
'
'P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).
'
'
'The arguments a and b are greater than zero, and x is
'nonnegative.
'
'ACCURACY:
'
'Tested at random points (a,b,x).
'
'               x     a,b                     Relative error:
'arithmetic  domain  domain     # trials      peak         rms
'   IEEE      0,1    0,100       100000      9.8e-15     1.7e-15
'   IEEE      1,5    0,100       100000      6.5e-15     3.5e-16
'   IEEE      0,1    1,10000     100000      2.2e-11     3.3e-12
'   IEEE      1,5    1,10000     100000      1.1e-11     1.7e-13
'
'Cephes Math Library Release 2.8:  June, 2000
'Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function FDistribution(ByVal a As Long, _
         ByVal b As Long, _
         ByVal x As Double) As Double
    Dim Result As Double
    Dim w As Double

    w = a*x
    w = w/(b+w)
    Result = IncompleteBeta(0.5*a, 0.5*b, w)

    FDistribution = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Complemented F distribution
'
'Returns the area from x to infinity under the F density
'function (also known as Snedcor's density or the
'variance ratio density).
'
'
'                     inf.
'                      -
'             1       | |  a-1      b-1
'1-P(x)  =  ------    |   t    (1-t)    dt
'           B(a,b)  | |
'                    -
'                     x
'
'
'The incomplete beta integral is used, according to the
'formula
'
'P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).
'
'
'ACCURACY:
'
'Tested at random points (a,b,x) in the indicated intervals.
'               x     a,b                     Relative error:
'arithmetic  domain  domain     # trials      peak         rms
'   IEEE      0,1    1,100       100000      3.7e-14     5.9e-16
'   IEEE      1,5    1,100       100000      8.0e-15     1.6e-15
'   IEEE      0,1    1,10000     100000      1.8e-11     3.5e-13
'   IEEE      1,5    1,10000     100000      2.0e-11     3.0e-12
'
'Cephes Math Library Release 2.8:  June, 2000
'Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function FCDistribution(ByVal a As Long, _
         ByVal b As Long, _
         ByVal x As Double) As Double
    Dim Result As Double
    Dim w As Double

    w = b/(b+a*x)
    Result = IncompleteBeta(0.5*b, 0.5*a, w)

    FCDistribution = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Inverse of complemented F distribution
'
'Finds the F density argument x such that the integral
'from x to infinity of the F density is equal to the
'given probability p.
'
'This is accomplished using the inverse beta integral
'function and the relations
'
'     z = incbi( df2/2, df1/2, p )
'     x = df2 (1-z) / (df1 z).
'
'Note: the following relations hold for the inverse of
'the uncomplemented F distribution:
'
'     z = incbi( df1/2, df2/2, p )
'     x = df2 z / (df1 (1-z)).
'
'ACCURACY:
'
'Tested at random points (a,b,p).
'
'             a,b                     Relative error:
'arithmetic  domain     # trials      peak         rms
' For p between .001 and 1:
'   IEEE     1,100       100000      8.3e-15     4.7e-16
'   IEEE     1,10000     100000      2.1e-11     1.4e-13
' For p between 10^-6 and 10^-3:
'   IEEE     1,100        50000      1.3e-12     8.4e-15
'   IEEE     1,10000      50000      3.0e-12     4.8e-14
'
'Cephes Math Library Release 2.8:  June, 2000
'Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function InvFDistribution(ByVal a As Long, _
         ByVal b As Long, _
         ByVal y As Double) As Double
    Dim Result As Double
    Dim w As Double

    
    '
    ' Compute probability for x = 0.5
    '
    w = IncompleteBeta(0.5*b, 0.5*a, 0.5)
    
    '
    ' If that is greater than y, then the solution w < .5
    ' Otherwise, solve at 1-y to remove cancellation in (b - b*w)
    '
    If w>y Or y<0.001 then
        w = InvIncompleteBeta(0.5*b, 0.5*a, y)
        Result = (b-b*w)/(a*w)
    Else
        w = InvIncompleteBeta(0.5*a, 0.5*b, 1.0-y)
        Result = b*w/(a*(1.0-w))
    End If

    InvFDistribution = Result
End Function


